﻿using System;

using FuncLib.DualFunctions;
using FuncLib.Functions;
using FuncLib.Optimization;
using FuncLib.Optimization.Ipopt;

namespace FuncLibSamples.Testing.Dual
{
	public class Graham
	{
		public static IPoint SampleSparse(int dimensions, bool useVectorAdd = false, bool flagVariables = false, bool indexedVariables = false, bool logging = false)
		{
			//Variable.VariableFactory vf = new Variable.VariableFactory();
			DualVariable[] v = new DualVariable[dimensions];
			VariableAssignment[] iv = new VariableAssignment[dimensions];
			double[] c = new double[dimensions];
			Function[] s = new Function[dimensions];
			Function[] r = new Function[dimensions];
			Function[] p = new Function[dimensions];

			Function sumS = 0.0;
			Function sumR = 0.0;
			Function sumP = 0.0;
			for (int i = 0; i < dimensions; ++i)
			{
				v[i] = new DualVariable("x" + i); // vf.CreateVariable();
				c[i] = 0.5;
				iv[i] = new VariableAssignment(v[i], 1.0);
				s[i] = 100 * Function.Pow(v[i], -2);
				r[i] = s[i] * v[i];
				p[i] = s[i] * (v[i] - c[i]);

				if (!useVectorAdd)
				{
					sumS += s[i];
					sumR += r[i];
					sumP += p[i];
				}
				if (flagVariables)
				{
					//s[i].Add(v[i]);
					//r[i].Add(v[i]);
					//p[i].Add(v[i]);
				}
			}

			/*if (useVectorAdd)
			{
				if (indexedVariables)
				{
					sumS = new Function.VectorAddFunction_IndexedVariables(s, dimensions);
					sumR = new Function.VectorAddFunction_IndexedVariables(r, dimensions);
					sumP = new Function.VectorAddFunction_IndexedVariables(p, dimensions);
				}
				else if (flagVariables)
				{
					sumS = new Function.VectorAddFunction_UsedVariables(s);
					sumR = new Function.VectorAddFunction_UsedVariables(r);
					sumP = new Function.VectorAddFunction_UsedVariables(p);
				}
				else
				{
					sumS = new Function.VectorAddFunction(s);
					sumR = new Function.VectorAddFunction(r);
					sumP = new Function.VectorAddFunction(p);
				}
			}*/

			double startR = sumR.Value(iv);

			// Prepare the optimizer.
			IpoptOptimizer o = new IpoptOptimizer();
			for (int i = 0; i < dimensions; i++)
			{
				o.Variables.Add(v[i]);
			}
			o.ObjectiveFunction = -sumP;

			for (int i = 0; i < dimensions; ++i)
			{
				o.Constraints.Add(0.5 <= v[i], v[i] <= 2.0);
			}

			for (int i = 1; i < dimensions; i += 2)
			{
				//Function f1 = (v[i] - v[i - 1]);
				Function f1 = (DualFunction.Disable(v[i]) - DualFunction.Disable(v[i - 1]));
				o.Constraints.Add(f1 == 0.2);
			}

			//o.Constraints.Add(sumR >= startR);

			if (logging)
			{
				// Verbose mode. Show Ipopt convergence.
				o.PrintLevel = 5;
			}

			IOptimizerResult or = o.Run(iv);

			Console.WriteLine(or.OptimalValue);
			return or.OptimalPoint;
		}    
	}
}
